
function DV=DVar(Z,s,ab,lambda,a,xi,n,w,nm)

N=length(s);

if w==1
    R=ab*exp(-lambda*Z);
else
    R=a*exp(-xi*Z);
end

S=diag(s);


if n>0
    DS=diag([zeros(n-1,1);1;zeros(N-n,1)]);
    DV=S*R*DS+DS*R*S;
elseif n==0
    DR=-Z.*R;
    DV=S*DR*S;
elseif n<0
    if w==1
        DR=R/ab;
    else
        DR=R/a;
    end
    DV=S*DR*S;
end

if ~isempty(nm)
    DV=DV(nm,nm);
end

